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AN OPTIMIZATION APPROACH FOR WELL-TARGETED 
TRANSCRANIAL DIRECT CURRENT STIMULATION 

SVEN WAGNER t, MARTIN BURGER *, AND GARSTEN H. WOLTERS t 


Abstract. Transcranial direct current stimulation is a non-invasive brain stimulation technique 
which modifies neural excitability by providing weak currents through scalp electrodes. The aim 
of this study is to introduce and analyze a novel optimization method for safe and well-targeted 
multi-array tDCS. For optimization, we consider an optimal control problem for a Laplace equation 
with Neumann boundary conditions with control and point-wise gradient state constraints. We 
prove well-posedness results for the proposed methods and provide computer simulation results in 
a highly realistic six-compartment geometry-adapted hexahedral head model. For discretization of 
the proposed minimization problem the finite element method is employed and the existence of at 
least one minimizer to the discretized optimization problem is shown. For numerical solution of the 
corresponding discretized problem we employ the alternating direction method of multipliers and 
comprehensively examine the cortical current flow field with regard to focality, target intensity and 
orientation. The numerical results reveal that the optimized current flow fields show significantly 
higher focality and, in most cases, higher directional agreement to the target vector in comparison 
to standard bipolar electrode montages. 

Key words, transcranial direct current stimulation, optimization, Neumann boundary control, 
gradient state constraint, sparsity optimal control, existence 

AMS subject classifications. 


1. Introduction. Transcranial direct current stimulation (tDCS) is a non-in- 
vasive, inexpensive and easy-to-perform brain stimulation technique which modihes 
neural excitability [23]. Changes in neural membrane potentials are induced in a 
polarity-dependent manner, for example, in a motor cortex study, anodal stimulation 
right over the motor cortex enhances cortical excitability, whereas cathodal stimula¬ 
tion inhibits it [23]. Recently, tDCS has been applied successfully in the treatment of 
neurological and neuropsychiatric disorders such as epilepsy [14], depression [6] and 
Alzheimer disease [13]. The effects of tDCS can be preserved for more than one hour 
after stimulation [24]. 

The conventional strategy is to apply the current density (TV<i> via two large elec¬ 
trodes, with the active electrode (anode) to be placed above the presumed target 
region and the reference electrode (cathode) far away from the target region [23, 24]. 
Accurate and detailed finite element (FE) head models have been created to inves¬ 
tigate the induced current density distribution [32, 10, 30]. While signihcant effects 
of stimulation as compared to sham were reported [6, 14, 23], computer simulation 
studies have revealed that the induced cortical current flow fields are very widespread 
with often strongest current density amplitudes in non-target brain regions [32, 30]. 
It is therefore a matter of debate whether the effects of stimulation are driven by the 
target brain region or elicited by adjacent cortical lobes [30]. 

In order to overcome the limitations of conventional bipolar electrode montages, 
algorithmic-based sensor optimization approaches were presented [10, 30, 21, 28]. Im 
and colleagues [21] searched for two electrode locations which generate maximal cur¬ 
rent flow towards a certain target direction. Ruffini and colleagues [28] described a 
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method for optimizing the configuration of multifocal tCS for stimulation of brain 
networks, represented by spatially extended cortical targets. Dmochowski and col¬ 
leagues [10] used a multi-channel array consisting of 64 fixed electrode locations to 
calculate optimized stimulation protocols for presumed target regions. They employed 
radial and tangential targets and reported that compared with conventional electrode 
montages their optimization approach achieved electric fields which exhibit simulta¬ 
neously greater focality and target intensity at cortical targets using the same total 
current applied. 

While the optimization results are very promising and a pilot study using multi-array 
tDCS devices for rehabilitation after stroke has been conducted [11], to our knowledge, 
there is currently no study providing an in-depth analysis of mathematical models for 
multi-array tDCS optimization. 

Given a volume conductor model D with a fixed electrode arrangement and a tar¬ 
get vector e in Dj with Dj C D being the target cortical area, the optimization 
approach estimates an optimal applied current pattern at the fixed electrodes. A 
Laplace equation with inhomogeneous Neumann boundary conditions is used to cal¬ 
culate the induced current density distribution [32] which is to be controlled by the 
boundary condition to ensure safe and focused stimulation. Therefore, the optimiza¬ 
tion problem for tDCS is in the class of control problems with Neumann boundary 
conditions [22]. 

In this paper, we modeled the electrodes with the point electrode model (PEM) 
[26, 1, 2]. They can, however, also be modeled using a complete electrode model 
(CEM) [26, 9, 12, 1, 2]. In [1, 2], it has been shown, that CEM and PEM only lead to 
small differences which are mainly situated locally around the electrodes and are very 
small in the brain region. Based on these results, the application of PEM is expected 
to result in negligible differences to the CEM and should thus provide a sufficiently 
accurate modeling of the current density within the brain region. 

This paper is organized as follows: In Section 2 we establish the mathematical model 
of tDCS and introduce the optimization problem for multi-channel tDCS. In Section 
3 we show existence of at least one minimizer to the optimization problem and the 
EE method is used for numerical discretization. In Section 4 we derive the algorithm 
for sensor optimization. Section 5 provides optimization results in a highly-realistic 
six-compartment head model (skin, skull compacta, skull spongiosa, CSF, gray and 
white matter) with white matter anisotropy. Furthermore, a conclusion and outlook 
section is presented. 

2. Mathematical Model of tDCS. In order to calculate the current flow field 
induced by tDCS, the quasistatic approximation to Maxwell’s equations is applied 
[32]. This yields the tDCS forward problem 


V • crV$ = 0 in D 
{aV^,n) = I onTcdQ 

$ = 0 ouTd = dn\T 

with $ being the electric potential, a G anisotropic conductivity 

tensor, I being the applied current pattern at the electrodes with non-zero values only 
at the electrode surfaces, n being the outward normal vector and P being a part of 
the boundary of the domain D. Furthermore, a Dirichlet boundary condition on the 
remaining part Pd is used to ensure that a solution to the tDCS forward problem is 
unique. We assume all parts of the boundary to have nonzero measure and to be of 
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Lipschitz regularity. We introduce the Sobolev space 

ij“5(r) :={uei?-5(r)| Ju{s)ds = 0 } c 

with being the standard Sobolev space for Neumann boundary values on F. 

The precise dehnition of H ~2 (F) is to be the dual space of (F) consisting of the 
Dirichlet traces of TF^-functions in ft. The integral in the above definition has to be 
interpreted as a duality product with the constant function equal to one, which is in 
iF 2 (F). The weak formulation of the forward problem is given by 

/ (tV$ • V'l' dx = / Id/ da, 

Jn Jr 

for all test functions G with vanishing Dirichlet trace on Fd. Note that the 

boundary integral on the right-hand side is only defined if / G T^(F), in general it is 
to be replaced by the duality product between FF “2 (F) and FFj (F). 

Under standard and naturally satisfied regularity assumptions, there exists a unique 
solution $ G FF^(fl) to the tDCS forward problem, which can be shown by the Lax- 
Milgram Lemma: 

Theorem 2.1. Let a G with a > agl^xs and let I G FF* ^(F). Then 

there exists a solution $ G H^{n) to the tDCS forward problem. The solution is 
unique ifTo has positive measure. Moreover, the following norm estimate holds: 


For current density optimization, a multi-channel array consisting of 74 fixed elec¬ 
trode locations (the locations of an extended 10-10 EEG electrode system) is used 
and optimized applied current patterns are calculated for presumed targets. Figure 
1 illustrates an overview of the optimization setup in a two-dimensional model. Op¬ 
timally, the induced current density distribution is maximal in the target region Lit 
and zero in non-target regions Llr = Ll\ Lit. Physically, a focal stimulation without 
stimulating non-target regions in not possible as the current has to flow through the 
volume to reach the target region. Therefore, the current density is restricted by e > 0 
such that 


|crV$| < e in 

Secondly, the total current applied to all electrodes is limited to 2 mA, a commonly 
used safety criterion [10]. 

For effective stimulation, the optimized current density distribution should be oriented 
perpendicularly to the cortical surface [5, 8]. Besides the correct target location, 
also the target direction is important as shown by Bindman and colleagues [5] and 
Creutzfeldt and colleagues [8] who were able to demonstrate in physical measurements 
in rats and cats, respectively, that the neural firing rate is strongly influenced by the 
direction of the field to the cortical surface. Perpendicularly-inwards (and parallel to 
the long apical dendrites of the large pyramidal cells in cortical layer V) stimulation, 
i.e. anodal stimulation, strongly enhanced the activity of the cortical neurons, whereas 
perpendicularly-outwards stimulation, i.e. cathodal stimulation, inhibited it. We thus 
maximize /^.^^(crVd*, e) da; with e C Lit being a unit vector perpendicularly-inwards- 
oriented to the cortical surface [5, 8]. 
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Fig. 1. Optimization setup in a two-dimensional model. The target region Q,t, a radial (sr) and 
tangential (et) target vector and the boundary F of the volume conductor model are demonstrated. 
The electrodes are depicted with solid lines. 


In the following we reformulate the current density constraint in the whole domain as 
wIctV^I < e with a weight w = 1 in O,. and 0 < a; ^ 1 in Qf Thus, we consider the 
constrained optimization problem 

(P) — / (ctV$, e) da; —min 

subject to w|crV<i)| < e 

J |-f|dx < 4 

V • crV$ = 0 infl 

(crV<i>, n) = I on r 
$ = 0 onTp) 

which is a control problem with Neumann boundary conditions [22]. Note that safety 
limitation on the inflow current is a control constraint, while the limit on the current 
density outside the target region is a state constraint. Also, the state constraint with 
w > 0 ensures that V$ is bounded in and due to the embedding theorem 

also bounded in This implies that <i> G and due to the trace theorem 

4)|r G Jl^(r). Since cr G jj^piigg that (o-V4),n) is bounded in i?“5(r). 

Remark 2.2. The constrained optimization problem (P) is a rather challenging 
optimal control problem. Firstly, the current density distribution tTV<l> is optimized, 
while many other applications require to optimize the potential field $. Seeondly, 
a point-wise gradient state constraint is used. While optimal control problems with 
state constraints were frequently used [17], not mueh attention was given to gradient 
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state constraints only until recently [29, 25, 33] focusing on distributed control in the 
domain however. Thirdly, we look for a control function I vanishing on large parts 
of the boundary, known as sparsity optimal control problem [19]. Indeed, we rather 
look for I as a combination of concentrated measures rather than an -function, so 
the above integral formulation has to be understood as formal notation for the total 
variation of the Radon measure identified with I. This will be made precise in the 
next section. 

Applying Gauss’ Theorem to (aV , e) dx provides a further motivation for the 
proposed objective functional. Assuming a and e locally constant we have 


/ (CTVd>,e)da;+ / tT$V • e + eVcr$ dx = / V • ((T$e) dx = / (a-<S>e,n}da 

J J J J dOt 

In order to maximize the current densities along the target direction, the potential 
should be maximized where the target vector e is oriented parallel to the outward 
normal vector n on dQ.t and minimized where it is antiparallel. 

3. Optimal Control Problem. The goal of this section is to provide insight 
into the existence theory of at least one minimizer to the tDCS optimization prob¬ 
lem (P). Furthermore, we show existence of at least one minimizer to a simplified 
optimization problem without control constraint at the electrodes. Finally, the finite 
element method is used for numerical discretization of the considered optimization 
procedure. 

3.1. Continuous formulation. We now reconsider the optimization problem 
(P). In order to provide a rigorous formulation of the safety constraint we interpret 
J as a Radon measure in the space At(F), whose norm, i.e. the total variation of a 
measure, is an appropriate replacement for the L^-norm. The constraint then becomes 


In order to obtain a unified formulation also allowing for other design constraints 

_ 1 

respectively a discrete electrode setup, we introduce a feasible set T’(F) C il<> ’^(F). 
To further simplify notation we define an operator A via 

A: ^ I^aV‘P. 

_ 1 

Since I G T>(r) C ^ (F), the boundary value problem for the Poisson equation has 
a unique solution $ G bounded by a multiple of the norm of I as guaranteed 

by Theorem 2.1. This implies that A is a well-defined linear operator. Hence we can 
substitute crV$ = AI in (P), which leads to the following equivalent reformulation 


(P?) — / (A/, e)dx—)• min 

Jn, levir) 

subject to w|AJ| < 6 

ll-^llTHCr) - 4- 

The existence of at least one minimizer I G T>(r) to (P^) is the topic of the following 
theorem. ^ 

Lemma 3.1. Let w be bounded away from zero in Q and let I G il* (F) be such 
that AI satisfies the constraint a;|A/| < <5 almost everywhere in LI. Then there exists 
a constant C depending only on the given data a, uj, S, LI and F, such that 


(3.1) 


'H~ 


qr) 


< c 
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Proof. We have 


By the trace theorem we can write each such as the Dirichlet trace of some 'I' G 
with vanishing trace on F^) and ||'I'||_f/i(n) < Ct for a universal constant in 
the trace theorem depending only on PL and F. Hence, by the weak formulation of the 
elliptic partial differential equation satisfied by $ with AI = (tV$, we find 

||r||^-i (P) < sup{ J^{AI) ■s/'b dx\ ll^'llffi(n) < Ct}. 

Using the Cauchy-Schwarz inequality and an elementary estimate of the L^-norm of 
AI in terms of its supremum norm, we find 


ll^ll 




<C = Ct 
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and we see that C only depends on ct, w, <5, U and F. □ 

Theorem 3.2. Let ViV) he weakly closed in ^ (F), let 0 G T>(F) and let e > 0. 
Then there exists at least one minimizer I G T’(F) of the convex optimization problem 

(P|)- 

Proof. An easy computation reveals that J = 0 G T’(F) is feasible for (P|), 

i.e. the set of feasible points is not empty. Moreover, it is closed in the weak-star 
_ 1 

topology of iFo ^ (F) n A^(F), which can be seen as follows: First of all the norm in 
Af(F) is lower semicontinuous with respect to weak-star convergence in this space, 
which implies that a limit of feasible points satisfies the safety constraint as well. 

Moreover, weak convergence of a sequence I in Hq ^ (F) implies weak convergence of 

AI^ in which implies that the limit satisfies the state constraint due to the 

weak closedness of pointwise constraints in Lf. Moreover, the feasible set is obviously 

_ 1 

bounded in M{Ll) due to the safety constraint and in ^(F) due to Lemma 3.1. 
Hence, the Banach-Alaoglu Theorem implies compactness in the weak-star topology. 
Since the objective is a bounded linear functional and in particular weak-star lower 
semicontinuous we obtain the existence of a minimizer by standard reasoning. □ 

We mention that our existence proof implicitely uses additional regularity induced 
by the gradient state constraint as done also in [33], however we do not discuss further 
regularity issues of the solution, which is more involved and possibly not to be expected 
with the nonsmooth terms in the objective functional. 

We observe that the objective as well as all terms involved in the constraints of 
the optimal control problem are one-homogeneous. This allows us to renormalize the 
unknown and in this way eliminate the safety constraint. This motivates the following 
simplified version: 


(Pg) — / (AJ, e)da:—>■ min 
Jn* /G'D(r) 

subject to wjAJj < e 

By renormalizing a minimizer of (Pg) we can obtain a minimizer for the full optimiza¬ 
tion problem (P^) with safety constraint for the control, which of course only holds 
if the total variation of the minimizer is finite: 
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Theorem 3.3. Let e > 0 be a threshold parameter and I G A^(r) n 77o ^(F) 


be a minimizer of the simplified minimization problem (Pe). Then I = 
minimizer of (P|) with 6 = 


41 


4€ 


l7M(r) 


ts a 




Proof. An easy calculation reveals that I = 


41 


fulfills the control and state 


constraint of (P^) with 5 = 


4e 


I At(r) 


N A 1 (r) 


. Assume that I is not a minimizer of (P|). Then 


there exists J* with smaller functional value, and —is feasible for (P|) with 
lower functional value than J, which contradicts the optimality of the latter. □ 

At a first glance the usefulness of Theorem 3.3 might appear limited, since the relation 
between the constraint bounds e and 5 is quite implicit, depending on the minimizer of 
the second problem. However, since there is no natural choice of the bound in either 
case one can choose 5 or e in appropriate ranges equally well. It turns out however that 
(Pg) is somewhat easier in particular with respect to numerical solutions. Moreover, 
it offers a better position to introduce sparsity constraints on the control by additional 
penalization, which is the subject of the following discussion. 


3.2. Penalized Problem Formulation. In order to control the applied current 
pattern at the fixed electrodes, the minimization problem (Pe) is extended by two 
penalties. While an L 2 term a da; is introduced to penalize the energy of the 
applied current pattern, an Li penalty f is used to minimize the number of 

active electrodes in the minimization procedure. The latter has to be reinterpreted in 
the same way as the safety constraint in terms of Radon measures. This leads to the 
following minimization problem 


{Pe’^) - [ {AI,e)dx+a [ /^da; + /3||7||^(p) ^ inin 

JUt Jr /GX>(r) 

subject to a;|AJ| < e 

We call the minimization problem (P“’^) with a > 0,/3 = 0 and a = 0,/3 > 0 to be 
the L 2 regularized optimization procedure (L2R) and the Li regularized optimization 
procedure (LIR), respectively. If both penalties are active, the penalty term is also 
known as elastic net. 

Since the penalty adds strict convexity to the problem, we can obtain uniqueness 
of a minimizer together with existence via analogous arguments as above: 

_ I, 

Theorem 3.4. Let 'DiT) be weakly closed in iJo ^, 0 G 'DiT), and let a > 0 
and 13 = 0. Then there exists a unique minimizer I G LfiT) to the L2R constrained 
optimization problem (P“’°). In the latter statement, /? = 0 can even be generalized 
to (3 >0. 

For the problem with pure L^-penalization the proof of existence is almost identical 
to the one of Theorem 3.3 and we hence conclude: 

Theorem 3.5. Let T) be weakly closed in ^ and let a = 0 and (3 > 0. 
Then there exists at least one minimizer I G AI(r) H UfT) to the LIR constrained 
optimization problem (P°Aj, 

3.3. Discrete Problem. For numerical discretization of the considered opti¬ 
mization problem we use the finite element method. The solution of the forward 
problem is approximated by looking for 4) = 'Yl!k=i with the being edge- 

based finite element basis functions vanishing on . 


(tV<I> • VT dx = / J 4' da, 
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for all in the span of {‘hfe}- Using standard techniques this allows to verify the 

_ 1 

existence and uniqueness of the discrete solution for arbitrary I G id* ^(r) and to 
introduce a discretized solution operator 

A : ido”®(r) ^ 


Let Tj, i = 1,..., JV denote the volume elements in the finite element discretization 
(tetrahedra or cuboids) and denote for a function u € the local mean value by 


b = TWTT / u(x) dx. 

\n Jt. 


Finally we discretize I = ^j4'j with (j)j being local basis functions on F. From 

the fact that I has mean zero we can eliminate Is- Now we introduce mappings Ai 
and A 2 incorporating the discretizations: 

s 

Ai:R^-\r)^V{r), Is = {Ii)^I = Y.h<ki 

A 2 : ((aV4>),)i=i,...,7v 


and define the discretized operator A 2 o A o Ai = B, represented by a matrix in 
]g3Arx(S-l)^ 

As the target vector e is only defined in the target region, we define a continuous 
mapping e: Qt 

e, ifT, cUt 
* \ 0, otherwise 


and replace {BIs,e) by {BIs,e), since {BIs,e) = 0 in U \ ilf We now introduce 
the discretized optimization problem (P“’^) as 

(P“’^) -{BIs,e)+a{Is,Is)+f3\\Ish^ min 
subject to uji\{BIs)i\ < e, 


where is a constant approximation of w in Ti (in our examples defined piecewise 
constant anyway). 

The existence (and potential uniqueness) of minimizers can now be verified exactly as 
for the continuous case above, defining an appropriate discrete version of the H~^- 
norm on 


s 

V{T) = {I = 

2=1 


e.g. via 

||/|| = sup{ J I<P da \ e span {$fc}, ||«'||//i(n) < C't}. 

For the sake of brevity we do not discuss the issue of, which can, e.g., be shown via 
F-convergence of the functionals involved. 
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4. Numerical Optimization. For numerical solution of the corresponding dis¬ 
cretized problem we employ the alternating direction method of multipliers (ADMM) 
[7]. The ADMM is a variant of the Augmented Lagrangian method, a class for solv¬ 
ing constrained optimization problems. The method combines important convergence 
properties (no strict convergence or finiteness is required) and the decomposability of 
the dual ascent method [7]. 

To solve the discretized optimization problem we substitute Big = y € 

and Is = z G and obtain the following Lagrangian y,z,Pi,P 2 ) which 

is to be minimized 

Li,^^f,^{Is,y,z,p^,P 2 ) = a{z,z) +I3\\z\\^ + ^{z-Is,z- Is) + {z - Is,Pi) 

-{y,e) + Y^y ~ Bis,y- Bis) + (y - bIs,p2) 
subject to uji\y^\ < e 

with p-i, /r 2 S K and p^ G P 2 G being the augmented Lagrangian parameters 

and the dual variables, respectively [7]. 

Is-step: Rearranging and neglecting terms without Is, putting the integrals together 
and expanding leads to a quadratic minimization problem 

z y,i 

+ ^(-2(y^ 57^+1) + {Blf\BI%+^) + —{Blf\pt)) 

Z P2 

Differentiating with respect to 7^"''^ and setting the equation system to be zero results 
in 


{paid + /r2S‘"B)7^+i - [yiz^ - = 0 

^ 7^+1 = {pild + p2B*^B)-\piz'^ - p 1 + P2i?‘’'y" - 

y-step: Rearranging and neglecting terms and putting the integrals together results 
in 


- 57^+1 - -pI - -e , y^+^ - BI^^^ - - -e) 

2 P2 fJ-2 P2 P2 

subject to WilyJ < e 

which can be solved analytically as follows 


= 


(bi'^+^+^pI+^&P 


rfe+l 


B)d’ 


if|(i?7^+^ + ip§ + ie),|>;^ 


(R7s+^ -t ^p| + otherwise 


U2 


z-step: Rearranging the terms, leaving out the terms without z and putting the 
integrals together leads to 






)+/ 3 ||^ 


/c+11 


2 


fe-l-l rfc-hl T 

Z - Ig ,Z - Ig 




Idz^+^ - 


Ml rk + l 


1%+^ - 




p\ 


k+l\^ 


- 7 


fc-i-i 


Pi) 


7^+1 - ^ 

•\/2pi 
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Algorithm 1 Algorithm for the discretized minimization problem (P“’^) 

1 
2 
3 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


Input: B,e,^Xi ,^2 , Ct, UJ^ G, Iprevi 7°, pO, pI z^, y°, N, TOL 
k = 0 


r/c _ T 
■*■8 ■‘■prev 


while k < 3 or 

T — 

prev — S 

jfc+l _ f,, I Utr 

for i = 1, - ■ ■ , 3N do 


> TOL do 


g = {pild + P2B^^B)-\piZ^ -p\+ 


Ai2 

— p 

y-2 ^ 


if 1 (57^+1 + + 


e {31%+^ + 

Vi 

|(B/g+W 

else 



(R7^+i + - 

end if 



— pn + —e)i 
tJ-2 


end for 

yk-\-\ - 


= {Id Idaid) ^Id 


Pi '^ = Li{p\ + - z 

= M + Bif 

k = k + l 

end while 

A — _4e_ 

o - II II 


Blf^ - y'^+i) 


fe+i\ 


111 -rk+l I 1 k 

+ TS^Pl ' 




P) 


19; return 


llAi(r) 

rk+l 

■'s ; 


<5, k. 


rfc+l 

-'s 


Ai(r) 


The solution to this equation is given as 

= i,Tu + .ur^,T(^ir + 

Dual update: Finally, the dual variables are updated 

p\+^=^,,{p\ + lf^-z^+^) 
pI+^ =^Ji2{pl + BIf^-y'^+^) 

We mention that a complete convergence analysis can be carried out following the 
arguments in [7]. Algorithm 1 depicts an overview of the optimization steps for the 
numerical solution of the discretized problem (P“’^). Note that the Euclidean norm 
was used for the stopping criterion in Step 3. 

5. Results and Discussion. We generated a highly realistic geometry-adapted 
six-compartment (skin, skull compacta, skull spongiosa, CSF, gray and white matter) 
finite element head model from a Tl- and a T2-weighted magnetic resonance image 
(MRI). For the compartments skin, skull compacta, skull spongiosa, CSF and gray 
matter we employed conductivity values a = 0.43S m“^, 0.007S m“^, 0.025S m“^, 
1.79 S m“^ and 0.33 S m“^, respectively [18, 27, 32]. For the modeling of the white 
matter anisotropy, a diffusion tensor MRI was used and the effective medium approach 
was applied. The effective medium approach assumes a linear relationship between 
the effective electrical conductivity tensor and the effective water diffusion tensor in 
white matter [31], resulting in a mean conductivity value a = 0.14S m“^ for the 
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Table 1 

Tangential target vector e; Averaged (over 924 targets) scaled threshold value 5, the number of 
iterations k and the scaling factor for tangential target vectors e and different input values 

e. 


e 

0.0001 

0.0005 

0.001 

0.005 

0.01 

5 [Ato“^] 

0.00332 

0.0085 

0.0154 

0.0175 

0.0208 

Iterations k 

132 

53 

33 

203 

143 


120.15 

233.96 

352.279 

1312.2 

1919.93 


white matter compartment [32]. As the current density amplitudes in the skin and 
CSF compartments are much stronger when compared to the brain compartments [32], 
we only visualized the current densities in the brain to enable best visibility of the 
cortical current flow pattern. The software package SCIRun was used for visualization 

[9]. 

Using this highly realistic volume conductor model we comprehensively examine the 
optimized current densities in the brain with regard to focality, target intensity and 
orientation. Firstly, the maximal current density in non-target regions is investigated 
for the discretized optimization problem (Pf’^) and four target vectors are used to 
calculate optimized stimulation protocols. For all simulations, the target area Uj 
always consists of the elements corresponding to the target vectors and only the 
volume conductor elements in the brain are used for optimization, i.e., w = 1 only in 
the brain compartments and w <i; 1 in the target region and in the CSF, skin and 
skull compartments. 

In order to investigate the current flow field that is induced by an optimized standard 
bipolar electrode montage, we use a maximum two electrode (M2E) approach. This 
approach stimulates only the main positive (anode) and the main negative electrode 
(cathode) of the LIR optimized stimulation protocols with a total current of 1 mA. In 
order to quantify the optimized current flow fields, we calculate the current densities 
in the direction of the target vectors (CDt) and the percentage of current density that 
is oriented parallel to the target vector {PAR), as shown in columns 4 and 5 in Table 
3, respectively. 

5.1. Maximal current density in non-target regions. In this section, the 
LIR approximated discretized optimization problem (P^’^) with (3 = 0.001 is used to 
calculate optimized stimulation protocols for a set of 924 tangential (parallel to the 
inner skull surface) and radial (perpendicular to the inner skull surface) target vectors 
and Theorem 3.3 is applied to estimate the averaged maximal current density in non¬ 
target regions as shown in Tables 1 and 2. Furthermore, the number of iterations until 
a minimum was found and the scaling factor ||F||^(r) is depicted. As can be seen 
in the tables, for the threshold value e = 0.001, the safety value 5 ensures that the 
maximal current density in the brain compartment is not dangerous for the subject. 
Furthermore, the number of iterations k until a minimum was found is lowest for e 
= 0.001. For this reason, a threshold value of e = 0.001 is thus used in this study as 
it allows safe and well-targeted stimulation and fast and robust computation of the 
stimulation protocol. Note that the number of iterations until a minimum is found 
increases with decreasing mesh size. However, because the maximal resolution of 
state-of-the-art whole-head MRI sequences is 1mm, we did not investigate mesh sizes 
smaller than 1mm. 
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Table 2 

Radial target vector e: Averaged (over 924 targets) scaled threshold value 6, the number of 
iterations k and the scaling factor for radial target vectors e and different input values e. 


e 

0.0001 

0.0005 

0.001 

0.005 

0.01 

S 

0.0043 

0.0114 

0.015 

0.023 

0.025 

Iterations k 

134 

36 

35 

68 

118 

A^(r) 

92.03 

176.08 

267.06 

887.08 

1597.73 


5.2. Mainly tangential target vector. Figure 2 displays the optimized cur¬ 
rent densities (for better visibility, we used a different scaling for the different rows) 
and the corresponding stimulation protocols (same scaling) for a mainly tangential 
target vector as shown in Figures 2-Al and -A2. As can be seen in Figures 2-Bl and 2- 
Cl, the optimized current flow fields using the L2R and LIR approaches, respectively, 
show high focality and maximal current densities can be observed in the target region. 
While the L2R and LIR approaches lead to target current densities of 0.022 A m“^ 
and 0.038 A m“^ (Table 3, second column), current densities in non-target brain re¬ 
gions are rather weak (Figs.2-Bl and -Cl). Overall, due to the rather widespread 
applied current pattern at the fixed electrodes (Fig.2-B2), the L2R optimized current 
flow field (Fig.2-Bl) has smaller amplitude when compared to the LIR one (Fig.2- 
Cl). On the other hand, the LIR stimulation protocol (Fig.2-C2) shows high focality 
with mainly two active electrodes, while only very weak compensating currents are 
injected at the neighboring electrodes. The M2E approach results in higher target 
current densities of 0.071 A m“^ (Table 3, second column). However, using the M2E 
approach, relatively strong current densities can also be noted in non-target regions, 
especially in pyramidal tract and deeper white matter regions (Fig.2-Dl). 

With a PAR value of 86.4, 86.8 and 85.9 for L2R, LIR and M2E (Table 3, fifth 
column), the current densities of all three approaches are mainly oriented parallel 
to the target vector, leading to CDt values, i.e., current densities along the target 
direction, of 0.019, 0.033and 0.061 A m“^, resp. (Table 3, fourth column). While 
the CDt values between LIR and M2E are thus less than a factor of 2 different, the 
averaged current density amplitudes in non-target regions is about 5.3 times higher 
when using the M2E approach (Table 3, third column). The LIR optimized current 
flow field thus shows significantly higher focality and also slightly better parallelity to 
the target vector in comparison to the bipolar electrode montage M2E. However, if no 
multi-channel tDCS stimulator is available, the M2E approach provides an optimized 
bipolar electrode montage for a mainly tangential target vector. 

5.3. Mainly radial target vector. The mainly radially oriented target is 
shown in Figures 3-Al and -A2. Figures 3-Bl and 3-Cl depict the optimized cur¬ 
rent densities for the L2R and LIR approaches, resp.. As shown in Table 3 (Column 
2), with a value of 0.045 A m“^, the target current density for the LIR approach is 
more than a factor of 1.7 times stronger than with the L2R approach (0.026 A m“^). 
Non-target regions show only weak current densities (Figures 3-Bl and -Cl). With 
a value of 0.063 A m“^, the M2E approach yields the largest target intensity (Table 
3, second column). However, for this approach, the maximal current density in the 
brain does not occur in the target region, but on a neighboring gyrus in between the 
stimulating electrodes and in mainly tangential direction (Figure 3-Dl). Moreover, 
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Fig. 2. A mainly tangential target vector (Figs. -Al and -A3): Figures -B -C and -D present 
optimization results for the L2R, LIR and M2E approach, respectively. The optimized current 
densities in the brain compartment (different scaling of rows), in a zoomed region of interest and the 
corresponding stimulation protocols (same scaling) are shown in Figures -1, -2 and -3, respectively. 


the M2E optimized current density shows again an overall much lower locality than 
the L2R and LIR current flow fields. 

The L2R stimulation protocol (Figure 3-B2) consists of an anode above the target, 
surrounded by four cathodes and a ring of very weak positive currents at the sec¬ 
ond neighboring electrodes, a distribution which might be best described by a sinc- 
function. As can be seen in Figure 3-C2, the LIR stimulation protocol is mainly 
composed of a main positive current at the electrode above the target region and four 
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Table 3 

Quantification of optimized current density. The averaged current density in the target area 
(CDa, second column), the averaged current density in non-target regions (third column), the inner 
product of current density and target vector (CDt, fourth column) and the percentage of current 
density that is oriented parallel to the target vector (PAR, fifth column) is displayed for different 
target vectors and methods (first column). 


Target 


[Am 

/n\nt dx 

J^{BIs,e) dx 

[%] 

l^tl 


1^4 

tangential L2R 

0.022 

0.00144 

0.019 

86.4 

tangential LIR 

0.038 

0.00151 

0.033 

86.8 

tangential M2E 

0.071 

0.0080 

0.061 

85.9 

radial L2R 

0.026 

0.00064 

0.025 

96.1 

radial LIR 

0.045 

0.00071 

0.043 

95.5 

radial M2E 

0.063 

0.0074 

0.048 

76.2 

patch L2R 

0.025 

0.00147 

0.022 

88.0 

patch LIR 

0.037 

0.00157 

0.033 

89.1 

patch M2E 

0.071 

0.0080 

0.062 

87.3 

deep L2R 

0.015 

0.00345 

0.013 

86.7 

deep LIR 

0.019 

0.00249 

0.018 

94.7 

deep M2E 

0.052 

0.01477 

0.049 

94.2 


return currents applied to the surrounding electrodes. 

The L2R and LIR optimized current flow fields show high directional agreement with 
the target vector e {PAR value above 95%), with the L2R slightly outperforming 
the LIR approach (Table 3, fifth column). With a PAR value of only 76.2% the 
current flow field in the M2E model shows much less directional agreement with the 
target vector. This results in CDt values of 0.025, 0.043 and 0.048 A m“^ for the L2R, 
LIR and M2E approaches, resp.. In order to obtain higher PAR values for the M2E 
approach, i.e., higher current densities along the target direction, the distance be¬ 
tween anode and cathode might be enlarged, in line with Dmochowski and colleagues 
[10] who reported that the optimal bipolar electrode configuration for a radial target 
vector consists of an electrode placed directly over the target with a distant return 
electrode. Another interesting bipolar electrode arrangement for a radial target vector 
might consist of an anode over the target and a cathodal ring around the anode. 

5.4. Extended mainly tangential target vector region. In this section, an 
extended target area of 3 mm x 1 mm x 1 mm is used for current density optimization. 
The target area is centered around the location of the superficial target vector that 
was also used in Section 5.2 and the corresponding target vectors are selected to be 
mainly tangentially oriented (Figures 4-Al and 4-A2). As can be seen in Figures 4- 
B1 and 4-Cl, the optimized current flow fields of the L2R and LIR approaches show 
high focality. L2R and LIR yield an averaged target intensity of 0.025 A m“^ and 
0.037 A m“^, resp. (Table 3, second column). Current density is not restricted to but 
focused to the target area, especially for the L2R approach (Fig. 4-Bl). In comparison 
to the superficial tangential target vector of Section 5.2, the LIR optimized averaged 
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Fig. 3. A mainly radial target vector (Figs. -Al and -A2): Figures -B -C and -D present 
optimization results for the L2R, LIR and M2E approach, respectively. The optimized current 
densities in the brain compartment (different scaling of rows), in a zoomed region of interest and the 
corresponding stimulation protocols (same scaling) are shown in Figures -1, -2 and -3, respectively. 


current density in the target area is decreased by about 3 % (from 0.038 A to 
0.037 A m“^) and the optimized stimulation protocols are also very similar. Because 
the main two electrodes are taken from the LIR optimization, this implies that also the 
M2E stimulation protocol remains constant to the stimulation protocol from Section 
5.2. In this way, the M2E approach yields an averaged target intensity of 0.071 A m“^ 
(Table 3, second column). 

With averaged PAR values of 88.0, 89.1 and 87.3 (Table 3, fifth column), the current 
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Fig. 4. An extended target region consisting of 3 tangential target vectors (Figs. -A1 and -A2): 
Figures -B -C and -D present optimization results for the L2R, LIR and M2E approach, respectively. 
The optimized current densities in the brain compartment (different scaling of rows), in a zoomed 
region of interest and the corresponding stimulation protocols (same scaling) are shown in Figures 
-1, -2 and -3, respectively. 


densities are mainly oriented parallel to the target vectors with LIR performing best. 
The averaged current flow field intensity along the target direction is 0.022, 0.033 
and 0.062Am“^ for the L2R, LIR and M2E approaches, resp. (Table 3, fourth 
column). However, for M2E also non-target regions reach significant current densities, 
as clearly shown in Eigure 4-Dl. The averaged current density amplitudes in non¬ 
target regions is 0.00157 and 0.0080 A m“^ for the LIR and M2E approaches (Table 
3, third column). The LIR optimized current flow field thus shows a factor of 5.1 
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Fig. 5. A deeper target vector (Figs. -A1 and -A2): Figures -B -C and -D present optimization 
results for the L2R, LIR and M2E approach, respectively. The optimized current densities in the 
brain compartment (different scaling of rows), in a zoomed region of interest and the corresponding 
stimulation protocols (same scaling) are shown in Figures -1, -2 and -3, respectively. 


higher locality in comparison to bipolar electrode montage M2E. However, the M2E 
approach provides an optimized bipolar electrode montage for an extended target area 
of tangential target vectors. 

5.5. Deep and tangential target vector. In the last simulation scenario, we 
investigate optimization for a deeper and mainly tangentially oriented target vector 
as shown in Figures 5-Al and -A2. Figures 5-Bl and -Cl depict the optimized cur¬ 
rent density distributions when using L2R and LIR for optimization, resp.. For those 
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approaches, the target current densities are 0.015 and 0.019 A m“^, resp.. With a 
value of 0.052 A m“^, which is more than 2.7 times the LIR value, the largest target 
intensity is, however, achieved with the M2E approach (Table 3, second column). 
CDt values of 0.013, 0.018 and 0.049 A m“^ lead to PAR values of 86.7, 94.7 and 
94.2 for the L2R, LIR and M2E approaches (Table 3, fourth and fifth column). The 
target current densities are thus for all three approaches oriented mainly parallel to 
the target vector. 

The L2R and LIR stimulation protocols show high locality with mainly two active 
electrodes, while only weak compensating currents are injected at the neighboring 
electrodes (Figures 2-B2 and -C2). Similar to Section 5.2, the compensating currents 
are stronger when using the L2R optimization procedure, leading to weaker target 
brain current densities as compared to the LIR stimulation. In order to enable cur¬ 
rent density to penetrate into deeper brain regions, the distance between the two main 
stimulating electrodes is larger when compared to the superficial mainly tangential 
target vector from Section 5.2, i.e., the electrode above the target region is not used 
for stimulation, while a more distant electrode is used as anode. 

For all three approaches, strongest current density amplitudes in the brain compart¬ 
ment always occur at the CSF/brain boundary above the target region (Figures 5- 
B1,-C1 and -Dl). This is due to the fact that the potential field V$ satisfies the 
maximum principle for harmonic functions which states that a non-constant function 
always attains its maximum at the boundary of the domain [15, Theorem 14.1]. Nev¬ 
ertheless, in average over all non-target regions, with a value of 0.014 77 A m~^ for 
M2E, the L2R (0.00345 A m“^) and LIR (0.00249 A m“^) optimized current flow 
fields show a factor of about 4.3 and 5.9 times lower current densities, resp. (Table 
3, third column). Overall, L2R and LIR thus have a much higher locality, which 
can also easily be seen in Figures 5-Bl,-Cl and -Dl. However, if no multi-channel 
tDCS device is available, the M2E approach provides an optimized bipolar electrode 
arrangement for a deep and tangential target vector. 

The deep target region does not seem to be located in a very deep region of the brain. 
It gets obvious that the deeper the target vector is located, the higher the averaged 
maximal current density in the brain compartment, especially in more lateral brain 
regions. Due to the maximum principle it is thus not possible to target deep regions 
without stimulating more lateral brain areas. However, many important target re¬ 
gions, e.g. auditory, motor or visual cortex, are located rather laterally, so that it will 
be possible to target with significant field strength in many applications. Moreover, in 
many applications of brain stimulation it might also not matter, if non-target regions 
are also involved, because the experimental setup focuses on the target region, for 
example when examining the change in event-related potentials (ERP) in pre- and 
post- tCS stimulation ERP measurements. 

6. Conclusion and Outlook. A novel optimization approach for safe and well- 
targeted multi-channel transcranial direct current stimulation has been proposed. Ex¬ 
istence of at least one minimizer has been proven for the proposed optimization meth¬ 
ods. For discretization of the respective minimization problems the finite element 
method was employed and the existence of at least one minimizer to the discretized 
optimization problems have been shown. For numerical solution of the corresponding 
discretized problem we employed the alternating direction method of multipliers. A 
highly-realistic six-compartment head model with white matter anisotropy was gen¬ 
erated and optimized current density distributions were calculated and evaluated for 
a mainly tangential and a mainly radial target vector at superficial locations, an 
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extended target area and a deeper mainly tangential target vector. The numerical 
results revealed that, while all approaches fulfilled the patient safety constraint, the 
optimized current flow fields show significantly higher focality and, with the exception 
of the L2R for the deep target, higher directional agreement to the target vector in 
comparison to standard bipolar electrode montages. The higher directional agree¬ 
ment is especially distinct for the radial target vector. In all test cases, because of 
a more widespread distribution of injected and extracted surface currents, the L2R 
optimization procedure (Pf’°) led to relatively weak current densities in the brain 
compartment. The LIR optimized current density distribution along the target direc¬ 
tion was in all test cases stronger than the L2R one and might thus be able to induce 
more signihcant stimulation effects. The stimulation will thus enhance cortical ex¬ 
citability especially in the target regions, while it will as good as possible prevent too 
strong excitability changes in non-target regions. 

We were able to demonstrate that the M2E approach provides optimized bipolar 
electrode montages as long as the target is mainly tangentially oriented. For radial 
targets, the M2E approach was unsatisfactory, an optimal bipolar electrode configu¬ 
ration might then consist of a small electrode placed directly above the target region 
with a distant return electrode or a small electrode over the target encircled by a ring 
return electrode, as proposed in [10]. 

A further application for the optimization method is transcranial magnetic stimulation 
(TMS). TMS uses externally generated magnetic fields to induce electrical currents to 
the underlying brain tissue [16]. Because there is no safety limit for the total currents 
applied to the stimulating coils but a safety-threshold for painful muscle twitching 
[16], the constrained optimization problem for multi-coil TMS is given as 

(Ptms) — / (crVd), e) da;—min 

Jnt 


subject to w < Em 



V • (tV$ = -V • 

dA{x, t) 

^ dt 

in n 

((tV$, n) = —{a 

dA(x,t) . 
dt ’ ^ 

on r 

$ = 0 

on Tjy 


with A{x,t) being the time-dependent magnetic vector potential and Em = 450 V m ^ 
being the threshold for painful muscle twitching [16]. By designing the changes in the 
magnetic vector potential one can consider as the optimization variables, re¬ 

spectively some parameters on which it depends linearly. The existence of at least 
one minimizer to the constrained optimization problem for TMS directly follows with 
similar arguments as in Theorem 3.2. 

Because the optimization method can be applied for both brain stimulation modali¬ 
ties, a combined tDCS and TMS optimization might outperform single modality tDCS 
or TMS optimizations, similar to what was shown for electro- (EEG) and magnetoen¬ 
cephalography (MEG) [3, 4]. While tCS is able to stimulate a radially oriented target, 
TMS is mainly not (like MEG is hardly able to detect radial sources [3, 4]). Possible 
applications of combined tDCS and TMS multi-channel and multi-coil optimization 
might thus be an improved stimulation of target regions containing both radial and 
tangential orientations or of deeper target regions. In order to induce action poten¬ 
tials in deeper target regions, the induced current density should exceed the threshold 
for neuronal depolarization of 150 V m“^ [16]. On the other hand, the threshold for 
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painful muscle twitching of 450 V must be kept [16]. The combination of the op¬ 
timized tDCS and TMS current density fields might lead to higher current densities 
in the target and simultaneously reduced current density amplitudes in non-target 
regions. 

While a thorough mathematical analysis of our novel multi-array tDCS optimization 
method was derived and results for different target regions were presented, besides the 
first promising results presented in [20], we did not yet further compare our method 
to the existing approaches in the literature such as, e.g., [10, 30, 21, 28]. Such a 
comparison is one of our future research goals. 
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